% convert model data to natural units, i.e. energy in units of hopping (t),
% times in units of t/hbar, distance in units of the lattice constant.

% input information
fname_input = 'expt_hydro_fit_bootstrap_test_summary.txt';

hopping = 925; %Hz = t/h. therefore, t/hbar = 2*pi*t / h
% hopping = 970; %Hz, for direction of modulation.
hopping_unc = 10; %Hz

% construct output filename
[~, name, ~] = fileparts(fname_input);
fname_output = fullfile([sprintf('%s_natural_units', name), '.txt']);

mli = 9.9883e-27; %kg
hbar = 1.0546e-34; %J*s
latt_const = 1064e-9 / sqrt(2); %m

%Read data
fid = fopen(fname_input,'r');
titles = split(fgetl(fid),'\t');
fclose(fid);

dat_mat = dlmread(fname_input,'\t',1,0);

% hbar_over_m = hbar / mli /(2 * pi * hopping) / latt_const^2;
% hbar_over_m_unc = hbar_over_m * hopping_unc / hopping;

Ts = dat_mat(:,1);
[Ts,sorted_i] = sort(Ts);
TsUnc = dat_mat(sorted_i,2);
Dens =  dat_mat(sorted_i,4);
DensUnc = dat_mat(sorted_i,5);

% gamma: MHz -> (t/hbar)
gamma_real = dat_mat(sorted_i, 7);
gamma_real_unc = dat_mat(sorted_i, 8);
gamma_nat = gamma_real * 1e6 /(hopping * 2 * pi);
gamma_nat_unc = gamma_nat .* sqrt( (gamma_real_unc ./ gamma_real).^2 +...
                                    (hopping_unc ./ hopping).^2 );

% deff: MHz * a^2 -> ta^2/hbar
deff_real = dat_mat(sorted_i, 9);
deff_real_unc = dat_mat(sorted_i, 10);
deff_nat = deff_real * 1e6/ (hopping * 2 * pi);
deff_nat_unc = deff_nat .* sqrt( (deff_real_unc ./ deff_real).^2 +...
                                    (hopping_unc ./ hopping).^2 ); 

% d_m: MHz * a^2 -> ta^2/hbar
dm_real = dat_mat(sorted_i, 11);
dm_real_unc = dat_mat(sorted_i, 12);
dm_nat = dm_real * 1e6 /(hopping * 2 * pi);
dm_nat_unc = dm_nat .* sqrt( (dm_real_unc ./ dm_real).^2 +...
                             (hopping_unc ./ hopping).^2 );
                         
% viscosity: MHz * a^2 -> ta^2/hbar
visc_real = dat_mat(sorted_i, 13);
visc_real_unc = dat_mat(sorted_i, 14);
visc_nat = visc_real * 1e6 /(hopping * 2 * pi);
visc_nat_unc = visc_real .* sqrt( (visc_real ./ visc_real_unc).^2 +...
                             (hopping_unc ./ hopping).^2 );

%save results
titles = {'T(t)', 'Tunc(t)', '<n>', '<n>_unc',...
    'gamma(t/hbar)', 'gamma_unc(t/hbar)',...
    'diff_eff(t/hbar/a^2)', 'diff_eff_unc(t/hbar/a^2)',...
    'diff_m(t/hbar/a^2)', 'diff_m_unc(t/hbar/a^2)',...
    'visc(t/hbar/a^2)', 'visc_unc(t/hbar/a^2)'};
dat_mat = [Ts, TsUnc, Dens, DensUnc,...
    gamma_nat, gamma_nat_unc, deff_nat, deff_nat_unc,...
    dm_nat, dm_nat_unc, visc_nat, visc_nat_unc];

if ~exist(fname_output, 'file')
    save_data_file(dat_mat, titles, '\t', '', '', fname_output, 0);
else
    warning('File %s exists. Did not overwrite. Manually delete if desired.');
end
